/* libwatopt++ - Optimization Toolkit
 *
 * Copyright (c) 2007 Christopher Alexander Watford
 * <christopher.watford@gmail.com>
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to
 * deal in the Software without restriction, including without limitation the
 * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
 * sell copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
 * IN THE SOFTWARE.
 *
 * $Id: sa.cpp 2 2007-05-06 18:01:15Z christopher.watford $
 */

#include <cstdlib>
#include <cmath>
#include <vector>
#include <limits>
#include <algorithm>

#include <opt.hpp>
#include <sa.hpp>

#include "utils.hpp"

namespace watopt
{

   ///
   /// SimulatedAnnealingOpt
   ///

   SimulatedAnnealingOpt::SimulatedAnnealingOpt(
	   double startTemp, size_t maxSteps, IAnnealingSchedule* sched,
      bool localOpt
	   )
   : _m_Temp(startTemp), _m_MaxSteps(maxSteps), _m_Schedule(sched)
   , _m_Best(NULL), _m_LastSteps(0), _m_LastTemp(0.0), _m_LocalOpt(localOpt)
   {
   }

   SimulatedAnnealingOpt::~SimulatedAnnealingOpt(void)
   {
   }

   void
   SimulatedAnnealingOpt::Initialize(IFitnessImplementation* guess)
   {
	   _m_Best = guess;
      seed_random();
   }

   IFitnessImplementation*
   SimulatedAnnealingOpt::Solve(void)
   {
	   if(NULL == _m_Best)
		   return NULL;

      if(NULL != _m_Instrument)
         _m_Instrument->Start();

      double tt = _m_Temp, curFit = std::numeric_limits<double>::infinity(),
         bestFit = std::numeric_limits<double>::infinity();
	   IFitnessImplementation* cur = NULL;
	   bool accept = false, conv = false;

	   _m_LastSteps = 0;
	   _m_LastTemp = tt;
	   _m_GoodAccepts = 0;
	   _m_MesaAccepts = 0;
	   _m_BadAccepts = 0;

	   // From step (0 .. _m_MaxSteps) perturb and anneal the guess
	   for(size_t kk = 0; kk < _m_MaxSteps; kk++)
	   {
		   // for now we exit on 0.0 instead of waiting until the end
		   if(bestFit == 0.0)
			   return _m_Best;

		   // cur is now a purturbed version of _m_Best
		   cur = _m_Best->Perturb();
		   if(cur == NULL) {
			   //ERROR!
			   return NULL;
		   }

         // if we want to use local optimization after each perturbation
         // invoke it now
         if(_m_LocalOpt)
         {
            IFitnessImplementation *localOpt = cur->LocalOpt();
            
            // LocalOpt can return NULL in error or if the fitness impl
            // doesn't support it
            if(NULL != localOpt)
            {
               delete cur;
               cur = localOpt;
            }
         }

		   // get our fitnesses
		   curFit = cur->Fitness();
		   bestFit = _m_Best->Fitness();
		   accept = false;

		   if(NULL != _m_Instrument)
            _m_Instrument->BeforeSelect(kk, tt, bestFit);

		   // Keep all better ones!
		   if(curFit < bestFit)
		   {
			   accept = true;
			   _m_GoodAccepts++;
		   } else if(curFit == bestFit) {
			   // Large 'Mesa' regions will have a big impact on the
			   // ability of simmulated annealing to find better
			   // solutions.
   			
			   //XXX: try 50/50 approach to mesa problems
			   if(rand01() < 0.5) {
				   accept = true;
				   _m_MesaAccepts++;
			   }
		   } else {
			   // Conditionally accept bad guesses
			   // Use the Metropolis Algorithm:
			   // rand < e^(-(cur - best) / (steps * temp))

			   double Metropolis = std::exp((bestFit - curFit) / (kk * tt));
			   if(rand01() < Metropolis) {
				   accept = true;
				   _m_BadAccepts++;
			   }
		   }

		   // did we accept the new object?
		   if(accept)
		   {
			   delete _m_Best;
			   _m_Best = cur;
		   } else {
			   delete cur;
		   }

		   // statistics
		   _m_LastSteps = kk;
		   _m_LastTemp = tt;

		   if(NULL != _m_Instrument)
            _m_Instrument->AfterSelect(kk, tt, bestFit, curFit, accept);

		   // anneal the temperature
		   tt = std::max<double>(_m_Schedule->Step(kk, tt, _m_MaxSteps, _m_Temp),
                                 0.0);
	   }

      if(NULL != _m_Instrument)
         _m_Instrument->Stop(bestFit);

	   return _m_Best;
   }

   void
   SimulatedAnnealingOpt::GetSolutions(
      std::vector<IFitnessImplementation*>& solutions) const
   {
      solutions.push_back(_m_Best);
   }
}
